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Calculations  were  performed  in  3-D  models  of  the  upper  75  -  ] 00  km.  of  the 
crust  ar.d  mantle  beneath  the  NORSAR  array  (-Thomson  and  Gubbins,  1982),  a  model 
for  a  region  in  northern  California  (Zandt ,  1981) ,  and  a  model  generated  by 
random  :  erturbations  to  a  1-D  velocity  structure  (-McLaughlin  and  Anderson,  1985) 

In  both  the  NORSAR  and  Zandt  models,  azimuthal  variations  in  teleseismic  amplitude 
were  found  to  be  on  the  order  of  a  factor  of  2  and  variations  in  travel  time  were 
found  to  be  on  the  order  of  several  0.1's  of  a  second.  These  models  had  a  maximum 
cf  4  to  9 ■  velocity  fluctuation  over  scale  lengths  of  10  to  100  km. 

A  significant  result  obtained  with  the  NORSAR  and  California  models  was 
that  the  scale-lengths  and  intensities  of  perturbations  were  such  that  the 
amplitude  variations  were  nearly  independent  of  frequency,  and  hence  adequately 
predicted  by  simple  ray  theory.  This  result  has  important  consequences  for  the 
yield  estimation  of  underground  nuclear  explosions  by  measurements  of  classical 
body  wave  magnitudes,  ,  versus  broader  band  measurements  of  radiated  energy 
in  the  time  and  frequency  domain.  If  deep  seated,  broad  scale  length  (50  km. 
and  areater) ,  velocity  anomalies  of  2%  or  more  are  a  common  occurrence  in 
the  upper  mantle  of  the  earth,  they  will  act  to  focus  and  defocus  body  waves  over 
a  broad  frequency  band.  The  focussing  and  defocussing  caused  by  these  broad 
anomalies  will  be  indq-endent  of  frequency  and  will  thus  introduce  a  scatter  in 
broader  band  measures  of  radiated  energy  which  will  be  equivalent  to  that  seen 
in  the  narrow  band  m^  measurement.  Focussing  and  defocussing  by  structure  in 
the  source  region  will  also  affect  the  coda  of  P  waves  if  a  portion  of  this  coda 
is  generated  in  the  receiver  again.  These  effects  may  help  explain  why  broader 
band  and  integrated  coda  measures  of  body  wave  energy  often  do  not  exhibit 
significantly  less  scatter  than  classical  m^  measurements. 

■  The  results  obtained  with  a  random  model  show  that  a  model  having  a 
maximum  velocity  fluctuation  as  small  as  0.8  percent  is  capable  of  producing 
caustics  and  multipaths  at  teleseismic  range.  The  production  of  multipaths 
strongly  depends  on  the  anisotrop>y  of  the  distribution  of  scale  lengths,  i.e., 
the  ratio  of  characteristic  vertical  and  horizontal  scale  length.  The  multipaths 
of  the  random  model,  however,  occurred  over  too  small  an  area  and  were  too  closely 
s: aced  in  arrival  time  to  be  resolved  with  standard  seismograph  systems  operating 
in  the  0.01  to  4  Hz.  band. 
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1  Introduction 


^ . 


The  3-D  structure  beneath  source  and  receiver  can  act  to  focus  and  defocus  teleseisrruc 
body  waves.  In  the  receiver  region,  the  focussing  and  defocussing  may  account  for 
variations  in  over  a  200  km.  aperture  array  as  large  as  that  seen  over  an  array 
having  teleseismic  dimensions.  The  following  sections  summarize  results  of  forward 
modeling  experiments  designed  to  measure  the  amplitude  fluctuations  predicted  by 
known  3-D  structure  beneath  sources  and  receivers  obtained  by  block  3-D  inversion  of 
travel  times.  The  importance  of  such  models  is  that  they  can  always  be  obtained  for  a 
particular  test  site  given  known  source  locations  and  times  within  a  test  site, 
supplemented  by  constraints  on  local  crustal  structure.  Given  such  models  and  the 
location  of  an  event  within  a  test  site,  it  is  possible  to  calculate  a  magnitude  bias  factor, 
which  can  be  used  to  correct  for  the  focussing/defocussing  effects  of  the  structure. 

This  factor  would  vary  as  a  function  of  event  location  within  the  test  site  and  azimuth  of 
the  receiver  station.  Similarly,  corrections  for  the  effects  of  structure  beneath  receiver 
arrays  may  be  formulated. 

This  report  details  the  results  of  the  first  complete  year  of  research  on  problems  of 
focussing  and  defocussing.  Several  source  and  receiver  models  have  been  investigated 
in  addition  to  a  model  of  NTS  by  Minster  ef  al  (1981).  The  scale  lengths  of  this  NTS 
model  were  too  broad  and  the  intensity  of  velocity  fluctuation  was  too  week  to  produce 
any  significant  amplitude  fluctuations.  Next  a  sequence  of  models  having  stronger 
velocity  fluctuations  were  investigated  in  order  to  determine  the  resolution  needed  to 
predict  significant  focussing  defocussing  effects.  Thus  far,  these  include  a  model 
beneath  NORSAR,  a  model  of  the  crust  and  uppermost  mantle  beneath  northern 
California,  and  a  random  model  constructed  to  satisfy  the  characteristic  statistics  of 
magnitude  variations.  The  results  reported  for  the  NORSAR  model  were  calculated  by 
Robert  N'owack  while  he  was  post-doctoral  fellow  at  M.J.T..  supported  by  this  contract 
The  experiments  with  the  /.andt  and  random  models  have  been  submitted  for 


publication  in  GJRAS  (Cormier.  1986).  A  paper  on  the  NORSAR  results  is  in  preparation 
by  Nowack  and  Cormier. 


2  Receiver  structure:  the  NORSAR  model 


2.1  DESCRIPTION  OF  MODEL 


The  structure  beneath  NORSAR  has  been  studied  by  a  number  of  investigators 
using  teleseismic  travel  time  data.  In  the  examples  shown  here,  the  model  derived  by 
Thomson  and  Gubbins  (1982)  is  used.  This  velocity  model  is  a  variation  of  the  model 
originally  derived  by  Aki  et  at  (1976),  which  was  the  first  application  of  seismic  travel 
time  inversion  for  a  3-D  structure. 


Good  quality  amplitude  data  has  been  recorded  at  NORSAR  and  has  been  studied  by 
several  investigators.  Kaddon  and  Husebye  (1978)  used  a  thin  lens  model  at  a  depth  of 
150  to  200  km  to  describe  both  the  amplitude  and  travel  time  data.  Statistical  models 
for  amplitude  fluctuations  at  NORSAR  were  investigated  by  Berteussen  (1975)  and 
Berteussen  et  at  (1975).  Thomson  and  Gubbins  (1982)  compared  the  amplitude  data  to 
that  predicted  by  the  travel  time  models  and  found  only  moderately  good  agreement, 
with  the  predictions  not  showing  large  enough  variations  across  the  array.  Thomson 
(1983)  attempted  a  separate  inversion  of  the  amplitude  data  with  again  only  a 
moderately  good  agreement  with  the  travel  time  inversion  results.  Thomson  (1983) 
suggested  the  possibility  that  ray  theory  may  not  give  reasonable  predictions  of  the 
amplitude  data  for  frequencies  of  1.0  to  3.0  Hz,  for  a  structure  such  as  NORSAR. 

Recently  an  amplitude  comparison  was  done  by  Haines  and  Thomson  (1986) 
between  the  Phase  Front  method  developed  by  Haines  (1983)  and  the  ray  amplitudes 
computed  using  the  ray  bending  method.  Their  comparison  showed  quite  different 
results  between  the  two  methods  and  their  interpretation  was  that  ray  theory  did  allow 
for  the  finite  frequency  effects.  This  is  a  very  interesting  result  since  the  NORSAR 
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velocity  model  is  smoothly  varying  with  length  scales  on  the  order  of  20.  km,  which 
approach  the  wavelengths  of  short  period  body  waves.  We  were  then  interested  in 
comparing  amplitudes  derived  from  ray  theory  with  Gaussian  beam  amplitudes  using 
this  NORSAR  velocity  model  C.  Thomson  provided  us  with  velocity  model  A2  from 
Thomson  and  Gubbins  (1982)  for  this  comparison. 

Figure  1  shows  this  velocity  model  for  NORSAR  The  model  is  180.  km  in  depth.  The 
spheres  in  Figure  1  represent  lower  velocities  and  the  cubes  represent  higher 
velocities.  Only  velocity  fluctuations  greater  than  3.5%  are  shown,  with  the  size  of  the 
symbol  corresponding  to  the  magnitude  of  the  velocity  fluctuation  from  the  average 
vertical  velocity  structure.  The  east  side  of  the  model,  with  values  of  x  greater  than  45. 
km,  is  in  general  fast  and  the  west  side  of  the  model  is  slow.  The  largest  fluctuations 
are  r  8%  and  occur  in  the  deepest  layer  from  96.  to  120.  km 

In  order  to  display  the  velocity  model  more  clearly.  Figure  2  shows  the  variations 
greater  than  1.5%  for  each  layer.  Layer  1  is  similar  to  Layer  2  with  smaller  relative 
velocity  fluctuations  and  is  not  shown.  Layer  2  is  from  24.  to  48.  km  and  has  higher 
velocities  in  the  south-central  regions.  The  Oslo  Graben  intersects  this  region  from  the 
south  and  trends  in  a  NNE  direction.  The  NORSAR  Array  occupies  the  middle  100.  km  on 
the  surface  of  the  model  region.  Layer  3  is  from  48.  to  72.  km  and  the  higher  velocities 
have  now  shifted  to  the  east.  Layer  4  is  from  72.  to  96.  km  The  velocity  fluctuations 
are  larger  with  lower  velocities  in  the  south  and  west  and  higher  velocities  to  the  east 
Finally,  the  deepest  layer,  Layer  5,  goes  from  96.  to  120.  km  the  inferred  depth  of  the 
lithosphere.  The  largest  fluctuations  are  seen  to  occur  in  this  layer  with  high  velociti'-- 
to  the  east  and  a  U-shaped  region  of  lower  velocities  to  the  south  and  west.  It  is 
interesting  to  consider  why  the  derived  velocity  fluctuations  increase  with  depth,  since 
the  resolution  of  the  derived  solution  of  these  lower  layers  appears  to  be  adequate  0n< 
might  intuitively  think  that  the  shallow  layers  should  be  more  heterogeneous  and 
therefore  have  a  greater  intensity  of  velocity  fluctuation.  Since  significant  velocity 
perturbations  continue  into  the  deepest  layer  of  this  model,  it  suggests  that  less 
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Figure  2.  Velocity  model  of  NORSAR  Variations  greater  than  1.5%  are  shown. 


intense  perturbations  that  continue  below  the  model  may  have  been  mapped  by  the 
inversion  into  the  lowest  layers  of  the  model. 

2.2  COMPARISON  WITH  THE  PHASE  FRONT  METHOD  AND  RAY  THEORY 

Figure  3  shows  the  location  of  the  NQRSAR  Array  and  three  specific  events,  with 
distances  from57=  to  70°.  used  by  Haines  and  Thomson  (1986).  To  compare  with  their 
results,  ray  amplitudes  using  the  shooting  method  were  computed.  To  do  this,  rapid 
two-point  ray  tracing  was  required  and  here  the  paraxial  ray  equations  were  used  to 
find  the  boundary  value  rays  (Cerven^  et  al  .  1984). 

Figure  4  shows  an  example  of  this  where  the  ray  from  A  to  B  is  given  and  we  want 
to  use  the  paraxial  ray  equations  to  determine  the  ray  from  A [  to  Bx.  Since  this  is  only 
an  approximation,  iteration  is  required  to  find  the  exact  ray.  The  problem  investigated 
here  is  a  slight  variation  of  this  where  an  initial  wavefront  is  given  with  some  initial  ray 
directions,  and  the  ray  eminating  from  the  wavefront  surface  going  to  a  particular 
station  is  required.  Using  the  iterative  application  of  the  paraxial  ray  equations, 
excellent  convergence  properties  were  found.  Typically  only  two  or  three  iterations 
were  required  to  get  within  .1  km  of  the  station  using  the  NORSAR  velocity  model. 

Figure  5  shows  the  rays  going  through  the  NORSAR  model  from  teteseismic  events 
A,  B,  and  C  in  Figure  3.  Twenty-five  equally  spaced  stations  were  specified  on  the 
surface  covering  the  general  location  of  the  NORSAR  array.  The  ray  trajectories  were 
checked  at  each  point  along  the  ray  by  using  the  eikonal  equation.  The  paraxial  ray 
equations  were  then  used  once  again  to  compute  the  amplitudes  at  the  stations.  Our 
derived  amplitudes  were  checked  with  ray  differencing  calculations  and  the 
comparisons  were  within  1.  percent. 

Our  computed  ray  amplitudes  were  then  compared  with  thfl  results  of  Haines  and 
Thomson  (1986).  In  comparing  the  ray  bending  amplitude  calculations  to  our  ray 
calculations,  there  were  some  discrepancies  at  particular  points  Part  of  this 


Application  of  paraxial  ray  equations 


Given  the  ray  Q(A,B),  find  the  approximate  ray  Q(A1,Bl). 

Application:  the  calculation  of  a  normal  ray 
from  an  initial  wavefront 


4  ^arax’aI  approximations  can  be  used  to  develop  an  iterative 
scheme  to  solve  the  two  point  ray  tracing  problem.  These  approxime 
tions  are  valid  in  the  neighborhood  of  ray  AB  and  can  be  used  to  find 
ray  Aj  Hj  that  solves  a  two  point  ray  tracing  problem. 


PARAXIAL  BOUNDARY  VALUE  RAY  TRACING 
EVENT  A 


PARAXIAL  BOUNDARY  VALUE  RAY  TRACING 
EVENT  B 


PARAXIAL  BOUNDARY  VALUE  RAY  TRACING 
EVENT  C 


Figure  5.  The  rays  from  the  NORSAR  model  from  the  teleseismic  events  A,  B, 
and  C  of  Figure  3. 


discrepancy  involved  boundary  value  rays,  shown  in  Figure  5,  that  exit  the  model  on  the 
sides  of  the  box.  Our  ray  amplitudes  resulted  from  extending  the  box  by  100.  km  on  all 
sides.  !t  is  unclear  how  the  ray  bending  calculations  handled  this.  These  particular 
events  had  distances  between  57c  to  70=.  For  the  actual  data,  distances  of  the  events 
ranged  from  25:  to  over  100’,  and  for  the  closer  events  this  problem  of  rays  exiting 
from  the  sides  would  increase. 

Excluding  the  stations  which  had  rays  that  exited  the  model  from  the  sides,  only 
four  points  for  the  three  events  A,  B,  and  C  varied  by  more  than  15%  in  amplitude  from 
the  ray  bending  calculations  of  Kaines  and  Thomson  (1986).  The  remaining  differences 
are  assumed  to  be  from  slight  parameterization  differences  of  the  model.  With  these 
considerations,  the  Phase  Front  calculations  of  Haines  and  Thomson  (1986)  were  not  in 
agreement  with  either  the  ray  bending  calculations  or  our  paraxial  ray  calculations. 
Also  the  differences  of  the  Phase  Front  calculations  with  frequency  an  the  ray 
calculations  did  not  vary  in  a  systematic  fashion.  Since  both  the  Phase  Front  and 
Gaussian  beam  methods  involve  similar  parabolic  approximations,  the  expectation  was 
that  they  should  have  been  in  excellent  agreement.  Haines  and  Thomson  note  that  the 
Gaussian  beam  method  would  be  inaccurate  in  the  deep  shadow  of  a  caustic  compared 
to  the  Phase  Front  method  The  NORSAR  model,  however,  generated  no  caustics  over 
the  short  ray  paths  of  the  expanded  plane  wave.  With  these  considerations,  our 
preliminary  conclusion  is  that  there  may  be  some  errors  in  Haines  and  Thomson's 
application  of  the  phase  front  method. 

2.3  FREQUENCY  DEPENDENCE 

In  the  following,  a  comparison  was  done  between  ray  theoretical  amplitudes  and 
Gaussian  beam  amplitude  for  several  frequencies.  To  simplify  the  comparison,  a 
vertically  incident  plane  wave  from  below  the  NORSAR  model  was  used.  Figure  6  shows 
the  variation  of  travel  time  at.  the  surface  of  the  model  for  the  vertically  incident  plane 


Figure  8.  Variations  in  travel  time  and  amplitude  for  a  plane  wave  incident  on 
the  NORSAR  model. 
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wave  from  below.  The  time  scale  is  1.0  sec.  Comparing  with  the  velocity  model  in  Figure 
2,  the  delayed  times  of  the  travel  time  surface  are  associated  with  the  slower  U  shaped 
velocity  region  in  the  lowest  layer  of  the  model  to  the  south  and  west.  The  earlier  times 
to  the  southeast  are  associated  with  faster  velocities.  The  corresponding  amplitude 
surface  for  this  case  is  also  shown  in  Figure  6.  The  amplitude  values  plotted  in  20.  Log 
(A/  S0).  are  shown  with  a  vertical  scale  going  from  -6.0  to  +6.0  db.  The  amplitude 
surface  can  be  seen  to  be  much  rougher  than  the  travel  time  surface.  This  is 
consistent  since  the  log  amplitudes  are  proportional  to  integrated  along  the  path. 
The  higher  amplitudes  are  seen  to  be  associated  with  the  delayed  travel  times. 

Next,  the  smaller  central  region  of  the  model  covered  by  the  NORSAR  seismic  array 
was  investigated.  The  amplitudes,  in  db.  for  this  smaller  region  are  shown  in  Figure  7. 
The  vertical  scales  go  from  -5.  to  +5.  db.  Berteussen  (1975)  found  at  NORSAR  a  range  of 
amplitudes  of  14.  db.  The  variation  shown  here  is  about  .6  db  for  this  vertically  incident 
case.  The  paraxial  ray  amplitude  surface  shown  in  Figure  7  has  a  high  peak  to  the 
northwest,  a  central  high  amplitude  region  with  a  lower  amplitude  to  the  south,  and  a 
gradual  increase  in  amplitude  to  the  southeast. 

The  Gaussian  beam  results  were  computed  using  a  simple  Gabor  wavelet  with  a  y  = 

6  and  several  center  frequencies.  Optimal  or  critical  beam  widths  are  specified  at  the 
source  plane  {Gerveny.  1962.  1985a,  1985b.  and  Klimes,  1984)  The  peak  amplitudes  ire 
then  plotted.  The  8.  Hz  Gaussian  beam  result  is  shown  in  Figure  7,  and  all  the  features 
shown  in  the  ray  amplitude  diagrams  are  present  But,  the  two  high  amplitude  regions 
are  slightly  lower  in  amplitude  than  the  corresponding  ray  results.  The  4.  Hz  Gaussian 
beam  result  shown  in  Figure  7  is  very  similar  except  the  amplitude  peak  to  the 
northwest  is  slightly  lower  still.  Finally,  the  1  0  Hz  Gaussian  beam  case  is  shown  and 
now  the  amplitude  peak  to  the  northwest  is  rounded  and  broadened  and  the  whole 
amplitude  surface  has  been  smoothed  Thus,  there  appears  to  be  a  progression  from 
the  ray  amplitude  results  to  the  high  frequency  beam  results,  with  greater  similarity  for 
higher  frequencies  The  frequency  dependence  between  1  to  8  Hz  ,  however,  is  small 


Rgure  7  in  amplitude  for  a  vertically  incident  plane  wave  on  the 

NuRSAR  model.  Results  for  several  frequencies  are  shown. 
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This  is  consistent  with  the  scale  lengths  of  this  model,  which  are  broad  (20  km.) 
compared  to  the  wavelengths  at  these  frequencies  (greater  than  8  km).  Since  the 
wavelengths  for  the  1.,  4.  and  8.  Hz  cases  are  8.,  2.,  and  1.  km  with  approximate  Fresnel 
zone  radii  of  20.,  12..  and  8.  km  the  Gaussian  beam  solutions  should  give  reasonable 
results  for  the  NORSAR  model,  at  ieast  down  to  1.  Hz. 

In  summary,  the  beam  amplitude  results  behave  in  a  similar  fashion  to  the  ray 
amplitude  results  with  the  lower  frequency  beam  solutions  becoming  smoother.  This 
suggests  the  possibility  of  using  lower  frequency  amplitude  data  to  smooth  over 
unwanted  amplitude  variability  to  obtain  results  for  an  equivalent  smooth  median.  The 
paraxial  ray  equations  have  been  used  for  the  2-point  boundary  value  ray  calculation 
and  have  shown  very  rapid  convergence. 

3  Source  structure 

3. 1  THE  NORTHERN  CALIFORNIA  MODEL 

The  focussing  and  defocussing  of  teteseismic  body  waves  by  3-D  structure  in  the 
vicinity  of  the  source  have  been  investigated  with  two  different  models.  The  first  model 
is  one  obtained  by  Zandt  (1981)  for  central  California  using  a  block  inversion  of 
teleseismic  travel  tunes  by  the  method  of  Aki  et  a l  (1976). 

The  Zandt  (1981)  model  has  four  layers  from  0.0  to  90.0  km  in  depth.  The 
horizontal  block  size  is  10.0  km  in  the  top  layer  and  20.0  to  25.0  km  in  the  lower  layers 
Average  velocity  variations  are  between  4,0  to  8  0  percent  in  the  top  layer  and  2.0  to  J.t 
percent  in  the  lower  layers.  The  rms  velocity  variation  measured  over  grid  points, 
however,  is  generally  much  lower,  on  the  order  of  less  than  2  percent.  This  is  because 
the  largest  variation^  trike  place  over  relatively  broad  regions,  having  characteristic 
scale  lengths  of  between  50  to  100  km.  (Figure  9) 
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Figure  8.  P  velocity  contours  in  a  horizontal  plane  at  90  kin.  depth  in  the  3-D 
model  for  central  California  by  Zandt  (1981)  Also  shown  are  the  pro¬ 
jections  of  source  locations  used  in  propagation  experiments  at 
teleseismic  range 


Seismograms  were  synthesized  in  this  model  by  summation  of  Gaussian  beams.  An 


explosive  point  source  was  assumed  at  the  center  of  the  model  at  9.6  km.  depth.  The 
3-D  model  was  patched  into  a  1-D.  flattened  whole  earth  model  by  use  of  the  propagator 
matrix  of  dynamic  ray  tracing  (Cormier,  1966).  The  elements  of  the  propagator  matrix 
n  were  calculated  in  the  3-D  source  region  by  numerical  integration  of  the  kinematic 
and  dynamic  ray  tracing  equations.  Velocities  and  their  first  and  second  order  spatial 
derivatives  in  the  3-D  region  were  defined  by  the  coefficients  of  cubic  spline 
interpolators  between  grid  points.  The  3-D  region  was  patched  into  a  1-D,  radially 
symmetric  earth  model  at  75  km  depth.  The  1-D  earth  model  was  the  1  Hz.,  isotropic 
PRFM  of  Dziewonski  and  Anderson  (1981).  PREM  was  first  flattened  using  the 
transformations  described  by  Muller  (1971)  The  FI  elements  were  then  computed  using 


dX 


a  fast  algorithm  m  which  the  quantities  and  X  are  given  by  a  analytic  formulae 


summed  over  thick,  vertically  inhomogeneous  layers  (Cerveny  and  Jansky,  1983).  TI  was 
determined  at  the  ray  end  points  in  the  receiver  region  by  propagator  multiplication 
The  matrices  Pp.  needed  to  evaluate  the  weighting  factor  for  superposition  of 
Gaussian  beams  are  given  by  the  2x2  sub-matrices  of  11 ,  and  H1Z  respectively. 
Focussing, Mefocussing  ’’fleets  of  the  structure  were  calculated  at  teleseismic  range  for 
the  variations  in  azimuth  and  variations  in  lateral  source  location  shown  in  Figure  3. 


3  11  FJffect  of  varying  azimuth  at  constant  source  location 

Figure  9  shows  the  results  of  beam  summation  for  a  teleseismic  P  wave  from  an 
explosive  source  embedded  in  the  Zandt  (1931)  model  At  source  location  sO, 
seismograms  were  computed  at  70°  for  eight  different  azimuths.  In  each  column  of 
Figure  9,  the  amplitudes  predicted  in  different  pass  bands  are  shown.  The  broadband 
pulse  was  that  obtained  using  the  source-time  function  of  Madariaga  and  Papadimitriou 
( 1985)  in  a  frequem  y  band  between  0.03  to  4  Hz.  Amplitudes  are  scaled  to  the 
maximum  peak  to  peak  amplitude  observed  at  the  101°  azimuth.  Amplitude  variations 
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Synthetic  seismograms  constructed  by  superposition  of  Gaussian 
beams  for  stations  at  70°  and  variable.-  azimuths  from  an  explosive 


are  on  the  order  of  two  and  travel  time  variations  on  the  order  of  several  0.1's  of  a 


second.  The  travel  time  variations  are  consistent  with  the  focussing/defocussing 
effects  --  large  amplitudes  correlate  with  slow  travel  times  and  small  amplitudes 
correlate  with  fast  travel  times.  The  largest  amplitudes  correlate  with  azimuths  in 
which  the  beams  sample  a  low  velocity  anomaly  to  the  southeast  of  source  sO.  This 
anomaly  is  a  significant  feature  in  both  layers  3  and  4  of  the  Zandt  model,  persisting 
over  50  km  of  depth  in  the  model.  Zandt  interprets  this  feature  as  a  \W-SE  trend  of 
lithospheric  thining  associated  with  a  fault  zone  that  includes  the  Calaveras.  Rogers 
Creek,  Maacama,  and  Lake  Mountain  faults. 

3  1  2  Effect  of  varying  source  location  at  constant  azimuth 


For  a  fixed  azimuth,  and  variations  in  source  site  from  positions  s30-  to  s30+ 
(Figure  10),  amplitude  variations  are  small  This  reflects  the  smaller  differences  in 
structure  between  the  regions  sampled  by  the  beams  compared  to  those  in  the 
azimuthal  experiment.  The  velocity  anomalies  in  the  deeper  layers  are  broad  features 
having  scale  lengths  of  50  km  or  more.  The  anomalies  in  the  shallower,  crustal  layers 
have  smaller  scale  lengths,  but  the  crustal  layers  are  thin  compared  to  the  total 
thickness  of  the  model  and  the  beams  spend  much  longer  time  in  the  thick  layers  3  and 
4.  Thus  the  broad  scale  lengths  of  the  anomalies  in  layers  3  and  4  have  the  greatest 
influence  on  amplitudes.  This  is  consistent  with  the  large  variations  in  amplitude  shown 
in  the  azimuthal  experiment  (Figure  9)  as  well  as  with  the  smaller  variations  in 
amplitude  due  to  changes  in  receiver  location  over  a  line  having  a  length  roughly  equal 
to  the  scale  length;  of  the  broad  anomalies  (Figure  10). 


V: 


3  i  3  Frequency  independence  and  Us  implications  for  treaty  verification 

The  relative  amplitudes  in  Figures  9  and  10  are  nearly  independent  with  respect  to 
the  frequency  band  of  observation  This  result,  has  important  consequences  for  the 


Figure  10.  Synthetic  seismograms  at  70°  distance  and  fixe< 
sources  at.  9.6  km.  depth  and  epicentral  locations  shown 


yield  estimation  of  underground  nuclear  explosions  by  measurements  of  classical  body 
wave  magnitudes,  m^.  versus  broader  band  measurements  of  radiated  energy  in  the 
time  and  frequency  domain.  If  deep  seated,  broad  scale  length  (50km  and  greater), 
velocity  anomalies  of  2%  or  more  are  a  common  occurrence  in  the  upper  mantle  of  the 
earth,  they  will  act  to  focus  and  defocus  body  waves  over  a  broad  frequency  band.  The 
focussing  and  defocussing  caused  by  these  broad  anomalies  will  be  independent  of 
frequency  and  will  thus  introduce  a  scatter  in  broader  band  measures  of  radiated 
energy  which  will  be  equivalent  to  that  seen  in  the  narrow  band  rn^  measurement. 
Focussing  and  defocussing  by  structure  iri  the  source  region  will  also  affect  the  coda  of 
P  waves  if  a  portion  of  this  coda  is  generated  in  the  receiver  region.  These  effects  may 
help  explain  why  broader  band  and  integrated  coda  measures  of  body  wave  energy 
often  do  not  exhibit  any  less  scatter  than  classical  measurements.  The  broadband 
and  coda  magnitudes  that  exhibit  the  least  scatter  typically  have  0.15  to  0.2  standard 
deviation  in  units  of  logarithm  of  energy  flux  rate  over  source  or  receiver  arrays  having 
apertures  of  200  km  (Bullitt  and  Cormier,  1984).  This  corresponds  to  about  a  1.5  to  2 
variation  in  the  amplitudes  of  particle  velocity,  similar  to  that  seen  in  the  synthetic 
seismograms  of  Figures  9  and  10.  These  results  suggest  that  knowledge  of  the  broader 
scale  length  velocity  anomalies  beneath  source  and  receiver  sites  may  be  useful  in 
correcting  and  reducing  the  scatter  in  magnitude  estimates  and  hence  the  uncertainty 
in  yield  estimates  of  nuclear  tests 

The  frequency  independence  of  the  amplitudes  calculated  in  the  Zandt  model  is  a 
characteristic  of  ray- theoret ically  predicted  amplitudes.  It  suggests  that  the  Gaussian 
beam  synthesis  should  not  be  necessary  to  accurately  calculate  amplitude  variations 
due  to  this  receiver  structure.  There  are  no  caustics  m  this  example  and  there  is  no 
particular  advaritug--  to  u'-.ng  Gaussian  beam  superposition  over  asymptotic  ray  theory, 
other  t  ban  exploiting  the  paraxial  <  stimat  ion  of  travel  time  Additional  evidence  of  why 
beam  superposition  ,n  this  example  reproduces  simple  ray  theory  is  illustrated  in 
Figure  1  1  Hay  densities  fi'.gur"  i  1)  w:t  hin  the  vicinity  of  a  receiver  can  accurately 
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Figure  1  1  The  end  points  of  rays  at  the  surface  of  the  earth  at  70°  at  high  and 
low  amplitude  stations  in  the  experiments  shown  in  Figures  9  and  10 
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predict  the  amplitude  observed  at  that  receiver.  Amplitude  is  nearly  proportional  to  1 
over  the  square  root  of  beam  density,  consistent  with  the  calculation  of  the  geometric 
spreading  by  ray  tube  area. 

All  combinations  of  sources  and  receivers  produce  ray  densities  that  are  similar  in 
form  to  Figure  11,  i.e..  uniform  over  broad  regions  surrounding  each  receiver,  with  no 
evidence  of  multipathing.  This  result  is  identical  to  that  obtained  by  Cormier  and  Aki 
(1982).  They  found  that  multipaths  are  not  generated  even  when  the  intensity  of 
anomalies  in  the  Zandt  model  are  doubled.  Multipathing,  however,  depends  on  the  scale 
length  of  anomalies  as  well  as  intensity.  This  is  emphasized  by  the  results  obtained  with 
the  second  model  investigated. 

3.2  A  RANDOM  MODEL 

The  second  model  was  one  generated  by  perturbing  a  1-D  model  at  10  to  20  km 
grid  points  in  horizontal  planes  with  an  rms  velocity  fluctuation  of  0.8%  (McLaughlin  and 
Anderson,  1985).  Unlike  the  Zandt  model,  this  model  introduced  caustics  and 
multipaths  at  teleseisrruc  range.  This  made  it  essential  to  use  Gaussian  beams  rather 
than  asymptotic  ray  theory  (ART)  to  synthesize  seismograms  at  receivers  in  the  vicinity 
of  caustics.  ART  evaluates  the  superposition  integral  by  a  stationary  or  saddle  point 
approximation.  The  stationary  phases  occur  at  the  discrete  rays  that  solve  the  two 
point  ray  tracing  problem  between  source  and  receiver,  leading  to  amplitudes 

proportional  to  factor  ^=====jy-.  This  factor  approaches  infinity  near  the  caustic 

surfaces  defined  by  detQ^  =  0  .  The  superposition  integral  and  its  integrand,  however, 
remain  regular  at  caustics  for  a  generalized,  complex  M  matrix  (fervent  et  a!.,  1992; 
terveny,  1985a, b) 


3  2  1  Effect  of  varying  source  Location  at  constant  azimuth 


Figure  12  shows  contours  of  velocity  at  the  bottom  of  the  random  model  and  the 
location  of  sources  in  a  variable  source  experiment.  The  velocity  contours  have  been 


left  unlabled  and  are  shown  only  to  illustrate  the  dramatically  different  scale  lengths  of 
velocity  fluctuation  in  this  model  compared  to  the  Zandt  model.  At  a  constant  azimuth, 
the  amplitude  fluctuations  (Figure  13)  due  to  variations  in  source  location  are  both 
larger  and  occur  more  rapidly  than  those  seen  in  the  Zandt  model  (Figure  10).  This 
reflects  the  fact  that  the  smallest  scale  length  of  velocity  fluctuation  (10  km.)  roughly 
equals  the  spacing  of  source  points.  Greater  frequency  dependence  of  the  amplitudes 
is  also  seen.  This  is  due,  in  part,  to  the  presence  of  caustics  in  the  vicinity  of  the 
receivers  for  some  of  the  source-receiver  paths.  Figure  14  is  a  plot  of  ray  end  points, 
illustrating  the  development  of  one  of  these  caustics  in  the  vicinity  of  the  70°  station 
for  a  source  at  location  slO-.  A  triplicated  zone  of  end  points  can  be  seen,  which  is 
elongated  along  a  narrow  azimuthal  zone.  Rays  having  end  points  wtthin  this  zone  are 
found  to  have  a  one  unit  advance  in  their  KMAH  index  (Ziolkowski  and  Deschamps, 

1980),  indicating  that  these  rays  have  passed  through  a  caustic  once.  A  receiver 
located  within  this  zone  of  triplicated  end  points  is  likely  to  observe  some  phase 
distortion  in  its  waveform  because  some  of  the  beams  that  contribute  to  the 

superposition  integral  will  have  a  phase  shift.  This  phase  distortion  is  difficult  to 

observe  in  synthetics  calculated  for  a  profile  of  stations  r400-  to  r400+  shown  in  Figur" 
lb  The  phase  distortion  appears  as  small  change  in  the  rise  time  of  the  broadband 
pulse  at  station  rO.  The  small  negative  first  break,  best  visible  on  the  broadband  pulse 
at  rO,  is  not  due  to  t.he  phase  shifted  beams,  but  rather  due  to  the  unique  interferenc¬ 
e-fleet  s  of  this  particular  beam  pattern.  The  largest  effects  to  be  observed  on  waveform- 
might  be  expected  on  the  short  period  instrument  The  highest  frequency  band  wonid 
have  the  highest  per  cent  ont  ribution  from  beams  within  the  triplicated  region  at  rb 
be  'a  use  beams  more  distant  from  rO  would  suffer  a  stronger  exponential  decay  No 
substantial  modifications,  however,  are  v.i  ;i  m  the  short  period  waveform  of  station  r 
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Figure  14  Rav  end  DOtnts  near  a  station  at  70°  distance  from  source  slO-  in 
the  random  model.  The  locations  of  a  profile  of  stations  in  the  vicinity 
of  a  caustic  intersection  is  given  by  the  tabled  points  r400-  to  r400+ 
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Comparison  of  amplitudes  in  different  frequency  bands  in  Figure  15  now  shows 
substantial  frequency  dependent  effects.  Long  period  amplitudes  vary  only  about  half 
as  much  as  short  period  and  broad  band  amplitudes.  A  much  broader  area  of  beams 
contribute  to  the  long  period  response,  smoothing  over  the  effect  of  the  caustic  region 
near  station  rO 

3  22  The  features  that  generate  caustics  and  multipaths 

Exactly  what  feature  of  the  random  model  was  responsible  for  the  thin,  elongated 
caustic  intersection  shown  in  Figure  14  ?  Since  this  feature  is  elongated  along  a 
particular  azimuth,  the  lateral  location  of  the  structure  is  constrained  to  be  along  the 
ray  paths  that  leave  the  3-D  portion  of  the  model  at  this  azimuth.  The  range  of  vertical 
take-off  angles  along  this  azimuth  is  also  bounded  by  the  apparent  edges  of  the  caustic 
intersections  in  the  plot  of  ray  end  points.  The  calculation  of  the  KMAH  index  can  also 
be  used  to  identify  the  particular  rays  that  are  tangent  to  the  caustic  surface  at  depth 
By  either  of  these  methods,  the  rays  that  describe  the  caustic  can  be  identified  and 
their  trajectories  plotted  through  the  3-D  region  of  the  model.  When  this  is  done 
(Figure  16),  it  can  be  seen  that  the  structure  responsible  for  the  caustic  at  teleseisnruc 
distance  is  a  low  velocity  zone,  extended  in  the  vertical  direction.  The  reason  why  such 
structures  have  been  generated  in  this  random  model  is  that  the  grid  spacing  at  which 
velocities  were  assigned  was  much  larger  in  the  vertical  direction  (30  km.)  than  in  the 
horizontal  direction  (10  or  20  km).  Thus,  there  will  occasionally  be  regions  of  the 
model  where  negative  perturbations  strongly  correlate  between  adjacent  vertical  grid 
lines,  forming  a  vertically,  elongated  zone  of  low  velocities.  Similarly,  elongated  zones 
of  high  velocities  will  be  formed.  The  surprising  observation  seen  with  this  particular 
model  is  that  only  a  very  small  perturbation  of  velocity  (0.8%),  with  a  vertical  scale  of  JO 
to  60  km.  and  a  horizontal  scale  of  10  km.  can  generate  caustics  and  phase  advances 
at  teleseismic  distances.  These  caustics  tntersect  the  surface  of  the  earth  at 
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teleseismic  range,  but  the  areal  extent  of  these  intersections  are  too  small  to  be  visibly 
identified  except  for  some  rather  subtle  effects  in  the  waveforms  at  a  few  number  of 
stations  The  generation  of  these  caustics  depends  on  the  strength  of  velocity 
fluctuation  as  well  as  on  the  relation  between  the  characteristic  vertical  and  horizontal 
scale  lengths  of  the  3-D  model.  In  the  example  considered,  the  shortest  scale  length  in 
the  vertical  direction  exceeds  that  in  the  horizontal  direction,  a  situation  which  is 
probably  not  the  common  state  of  crust/lithospheric  structure  (McLaughlin,  personal 
communication).  Some  notable  exceptions  to  this  would  include  regions  having 
concentrations  of  intrusive  pipes  and  plumes.  The  results  of  the  single  modeling 
experiment  described  here  suggest  that  some  distributions  of  heterogeneity  would 
produce  unacceptably  large  effects  on  teleseismic  waveforms.  It  is  clear  that  forward 
modeling  of  the  effects  of  very  general  distributions  of  heterogeneities  will  be  useful  in 
defining  the  "heterospectrum''  of  the  lithosphere.  (The  heterospectrum  is  a  term 
adopted  by  Wu  [  1996]  to  embrace  both  the  the  strength  of  velocity  fluctuation  and  its 
three-dimensional  spatial  spectrum). 

4  Validity  of  the  results 

All  the  example  seismograms  were  calculated  by  superposition  of  Gaussian  beams 
The  accuracy  of  this  technique  depends  both  on  the  validity  of  using  the  first  term  in 
an  asymptotic  solution  to  the  elastodynamic  wave  equation  and  making  a  Taylor 
expansion  of  the  complex  phase  about  the  central  ray  (the  paraxial  approximation)  It 
is  thus  appropriate  to  question  whether  the  3-D  models  of  velocity  discussed  in  this 
paper  have  exceeded  the  domains  of  validity  of  these  approximations.  Validity 
constraints  in  continuous  media  have  been  formulated  by  Bcydoun  and  Ben-Menahem 
(1995)  and  Cerven^’  ( 1 995b)  White  et  al  (1996)  have  considered  problems  encountered 
with  continuous  media  as  well  with  boundary  interactions  The  simplest  constraints  »o 
check  are  those  related  to  the  asymptotic  expansion,  which  assumes  decoupling  of  P 


and  S  waves  to  zeroth  order  and  neglect  reflections  and  conversions  by  regions  of 
strong  gradient.  These  constraints  require  that  wavelength  be  much  less  than 


quantities  such  as  and  (Kravtsov  and  Orlov.  1980;  Beydoun  and  Ben-Menahem. 

1985).  where  v  and  p  are  velocity  and  density  respectively.  Both  the  Zandt  model  and 
the  random  model  satisfy  this  constraint  throughout  the  frequency  band  0.03  to  4  Hz. 


Constraints  on  the  validity  of  the  paraxial  approximation  and  the  beam 
superposition  have  been  cast  in  terms  of  the  real  and  imaginary  parts  of  the  II  matrix 
and  distance  from  the  central  ray  (Beydoun  and  Ben-Menahem.  1985;  Cerven^,  1985b). 
In  order  to  see  if  these  constraints  are  obeyed  everywhere  along  a  beam,  the  H  matrix 
must  be  known  everywhere  along  the  central  ray  associated  with  that  beam  At  any 
point  S  along  the  central  ray.  11(5)  can  be  found  by  using  the  II  matrix  to  back 
propagate  the  complex  M(Ct)  matrix  selected  at  the  end  point  O,.  Although  these 
constraints  were  not  calculated  in  the  two  examples  described,  a  check  on  the  overall 
validity  of  the  beam  superposition  was  taken  from  the  results  of  reciprocal 
experiments,  in  which  the  positions  of  sources  and  receivers  were  reversed.  Examples 
of  this  test  in  2-D  media  are  given  by  Nowack  and  Aki  (1984)  and  Muller  (1984). 
Reciprocal  experiments  were  conducted  in  the  form  of  allowing  a  plane  wave  to  be 
vertically  incident  on  the  3-D  models.  The  plane  wave  was  expanded  into  Gaussian 
beams  by  the  procedure  described  by  Cerven^  (1982).  and  the  wavefield  was  calculated 
at  receivers  on  the  surface  of  the  3-D  models.  Such  an  experiment  was  conducted  on 
the  random  model  (McLaughlin  and  Anderson,  1985)  and  on  a  model  having  a 
heterospectrum  similar  in  scale  length  and  velocity  fluctuation  to  the  Zandt  model 
(Nowack  and  Cormier,  1985)  Although  the  geometry  was  not  precisely  reciprocal  and 
waveforms  were  not  directly  compared,  both  experiments  produced  intensities  and 
coherence  of  amplitude  fluctuations  similar  to  those  produced  by  the  teleseismic 
experiments. 


5  Conclusions 


n . 


The  theory  and  examples  discussed  in  this  paper  have  shown  how  the  propagator 
matrix  of  the  dynamic  ray  tracing  equations,  TI ,  can  be  exploited  to  connect  3-D  to  1-D 
portions  of  a  model  Both  plane  wave  and  point  source  initial  conditions  are  required  to 
specify  the  elements  of  the  propagator  matrix.  Thus,  both  plane  wave  and  point  source 
solutions  are  generally  useful  for  all  asymptotic  methods  of  body  wave  synthesis  and 
not  just  for  Gaussian  beams 

An  investigation  was  made  of  the  effect  of  3-D  structure  in  the  source  region  on  the 
focussing  and  defocussing  of  teleseismic  P  waves.  These  effects  were  observed  with  two 
different  models  as  a  function  of  source  location  within  the  3-D  model  and  of  azimuth  at 
the  source.  In  a  3-D  model,  block  inverted  from  teleseismic  travel  times,  ray  theoretical 
amplitudes  matched  those  predicted  from  superposition  of  Gaussian  beams.  In  this 
model,  the  characteristic  scale  lengths  of  the  most  intense  velocity  fluctuations  (4”) 
were  on  the  order  of  50-100  km  in  the  source  region.  This  model  produced  a  factor  of 
two  fluctuation  in  amplitude,  associated  with  fluctuations  in  travel  time  on  the  order  of 
several  tenths  of  a  second.  Amplitude  variations  due  to  variations  in  source  location 
were  small  over  location  variations  that  were  small  with  respect  to  the  scale  length  of 
velocity  fluctuation  .All  amplitude  variations  were  nearly  independent  of  frequency 
across  the  body  wave  band  of  0.03  to  4  Hz.  The  frequency  independence  of  amplitude 
variations  across  the  body  wave  band  may  have  important  implications  for  removing  the 
effects  of  azimuthal  amplitude  variations  due  to  3-D  structure  beneath  nuclear  test 
sites  Broad  scale  length,  deep  seated  structure  can  affect  short  period  as  well  as  broad 
band  and  coda  measures  of  radiated  seismic  energy,  its  effects,  however,  may  be  easily 
correctable  if  a  3-D  model  of  the  source  region  is  known  from  block  inversion  of  travel 
time  residuals.  A  resolvable  block  size  of  about  20  km  may  be  all  that  is  necessary  to 
formulate  corrections  based  on  azimuth  of  teleseismic  station  and  source  location 


Calculations  with  a  random  velocity  model  demonstrated  that  a  smaller  intensity  of 
velocity  fluctuation  (0.8%)  can  produce  even  larger  amplitude  variations  if  the  smallest 
scale  length  of  fluctuation  is  on  the  order  of  10  km.  The  random  model  was  constructed 
such  that  the  smallest  scale  length  of  velocity  fluctuation  was  shorter  in  the  horizontal 
direction  than  in  the  vertical  direction.  At  teleseismic  range,  this  model  generated 
isolated  caustics,  which  were  elongated  along  the  azimuth  of  approach.  The  waveform 
distortion  associated  with  these  caustics  was  small.  The  fact  that  the  caustics  were 
generated  at  all  by  such  mild  3-D  perturbations  is  significant.  It  suggests  that  this  type 
of  synthetic  modeling  may  be  useful  in  limiting  some  of  the  attributes  of  the 
heterospectrum  of  the  earth's  lithosphere. 
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